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Abstract 

We study QCD thermodynamics using two flavors of dynamical overlap 
fermions with quark masses corresponding to a pion mass of 350 MeV. We 
determine several observables on N t = 6 and 8 lattices. All our runs are per- 
formed with fixed global topology. Our results are compared with staggered 
ones and a nice agreement is found. 
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1. Introduction 

At high temperatures the dominant degrees of freedom of strongly inter- 
acting matter change from hadrons to quarks and gluons. This transition 
can be studied using lattice gauge theory. There are various results using 
different fermion regularizations. 

Most of these results [H, 0, [si use the computationally least expensive 
staggered discretization which also preserves some of the continuum chiral 
symmetry. Even though different staggered results seem to agree with each 
other one should not forget that all these works use the fourth root trick 
to study Nf = 2+1 flavors of quarks. There is still an ongoing debate 
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in the literature about the correctness of this approach. Furthermore taste 
symmetry breaking may lead to large discretization errors when using small 
quark masses, especially at low temperatures. 

There are also several results using Wilson fermions 0,0, Ifij]. Since Wilson 
fermions break chiral symmetry explicitly one has to take very fine lattices 
to study chiral symmetry restoration at finite temperature. Due to the scat- 
tering of the low lying eigenvalues of the Wilson-Dirac operator one needs 
large lattice volumes when going to small pion masses. There are also first 
results with twisted mass fermions [7J. 

It seems logical to use chiral fermions to study chiral properties at finite 
temperature. Even though lattice chiral fermions are computationally much 
more expensive than the other types of discretization there are results in the 
literature using domain wall fermions Js|] as well as first attempts with overlap 
fermions j^] . While do main- wall fermions provide exact chiral symmetry onh 



for an infinite extent of the fifth dimension, the overlap formulation [10 



has the advantage of exact symmetry on finite four dimensional lattices |l2j . 

In this exploratory Letter we present results using two degenerate flavors 
of dynamical overlap fermions. We use two lattice resolutions corresponding 
to Nt = 6 and 8 temporal extents. We determine the temperature dependence 
of the chiral condensate, the chiral susceptibility, the quark number suscepti- 
bility and the Polyakov loop. The results are compared to Nf = 2 staggered 
data, the details of these simulations are summarized in Appendix B 



2. Lattice action and simulation details 

The possibility of using the Hybrid Monte Carlo algorithm (HMC) with 
overlap fermions was first discussed in Reference (l3| . The overlap operator 
was implemented with a multi-shift inverter using the Zolotarev rational ap- 
proximation [14| • The working of the algorithm was demonstrated on a small 
thermodynamic study on 6 3 ■ 4 lattices. It was also observed that treating 
topology changes requires special care during the HMC trajectories. One 
has to track the lowest lying eigenvalues of the Wilson kernel of the overlap 



operator. This was studied in detail in References [15|, |16|, ll7| • It was demon- 



strated in Reference [18j that one can do simulations with a fixed topological 
charge in several different sectors and it is possible to determine their relative 
weight. However, even this approach requires a tracking of Wilson eigenval- 



ues. In Reference [19( it was suggested that by adding an extra heavy Wilson 
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fermion to the action which decouples in the continuum limit, one can dis- 
able topological sector changes and at the same time speed up the algorithm 
significantly. It was also claimed that in the thermodynamic limit physics is 
independent of the global topology and therefore this approach should give 
correct results. However, significant power-like finite volume corrections are 



expected [20|,[21(. Here we follow the same approach: we add an extra Wilson 
fermion to suppress low lying eigenvalues of the Wilson kernel and disable 
tunneling between different topological sectors. As a further improvement 
we use smearing in the Wilson kernel. It was observed in jiij that smearing 
significantly improves the properties of the overlap operator. Furthermore, 
since smearing decreases the eigenvalue density in the middle of the Wilson 
spectrum it results in a significant speedup of the algorithm (ijf . 

In the gauge sector we use a tree level Symanzik improved gauge action. 
The overlap operator can be written as 

/ Tfl \ 

D= [m -—j (l + 7 5 sgn(# H /))+m, (1) 

where Hw = 75-Djy is the Hermitian Wilson operator with a negative —2 < 
—mo < mass parameter and m is the mass of the overlap quark. For the 



Wilson kernel we use two steps of HEX smearing [24], [25j, [26[ with smearing 
parameters of a>i = 0.72, «2 = 0.60 and a 3 = 0.44. In order to set m we 
evaluated the Wilson kernel on quenched configurations with the targeted 
lattice spacings in this work and located the point which is in the middle 
between the physical modes and the first doublers. This resulted in mo=1.3. 
The simulations are performed with Nf = 2 flavors. 



As suggested in 19] we add two irrelevant terms to the action to suppress 



low eigenvalues of H w and fix topology: 

Se = 2j {i'E{x)D w (-mQ)i) E (x) + (f>\x)[D w (-m ) + im B -f 5 T 3 }(p(x)} . 

X 

(2) 

The first term is the action of two flavors of extra fermions with negative 
mass — m . The second term, including a two component bosonic field is 
included to control the effect of the extra fermions. The eigenvalues of H\y 
below itlb are most strongly suppressed. Since both mo and uib are fixed in 
lattice units they correspond to infinitely large masses in the continuum limit 
and both these terms decouple. For the bosonic mass we use ms=0. 54. Since 
our lattice action results in a fixed topology we aimed to make simulations 
with zero topological charge. 
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Figure 1: Left: The lattice spacing as a function of the (3 coupling. The opened box shows 
a run with a fixed topological charge of Q = 20. Right: the bare quark mass as a function 
of p. 



We use a HMC algorithm with the Hasenbusch trick |27J , with an Omelyan 
integrator (28| and with a Sexton- Weingarten multi-scale scheme for the dif- 
ferent fields |29j . The latter ingredient turned out to be rather advantageous, 
the extra Wilson fermion has to be integrated with a much smaller stepsize 
than the much more expensive overlap fermion. 

The first step of our analysis was to determine the line of constant physics 
(LCP) and the scale. We used 12 3 • 24 lattices for = 3.6,3.7,3.8 and 3.9 
and 16 3 • 32 for j3 = 4.0 and 4.1 with some initial guesses of the bare quark 
masses between 0.02 and 0.06. We determined the wo scale [3o[ as well as the 
pion masses on all of these lattices. The wq scale is defined implicitly by the 
equation d/dt[t 2 E(t)]\ t=w 2 = 0.3, where E(t) is the expectation value of the 
gauge action evaluated on configurations evolved by the Wilson flow 3l| with 
parameter t. Since due to the chiral symmetry of overlap fermions m 2 oc m 
and the Wq scale is quite insensitive to the quark mass, it was possible to tune 
the quark masses to have a fixed value of m n •Wo=0.312 for each beta without 
further simulations. With the physical value of Wq at the Nf = 2 + 1 flavor 
physical point, w = 0.1755 fm, this would correspond to m n = 350 MeV. 
Rigorously a conversion into physical units is only well defined in QCD with 
physical quark masses, so the results in MeV or fm are for orientation only. 
The lattice spacing as a function of (3 and the LCP are shown in Figure [TJ 
In one of our runs we had a nonzero topological charge, Q = 20. One can 
see from the plot that both the scale and the LCP have still a significant 
dependence on topology for our volumes. 
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For the finite temperature runs we use two sets of lattices: 12 3 • 6 and 
16 3 • 8. Since the lattice spacing and the LCP is quite ambiguous at large 
lattice spacings and the algorithm performs poorly on coarse lattices we 
decided not to do N t = 4 runs. We had 9 different (3 values for both N t = 6 
and 8 in the range (3 = 3.6 . . . 4.1. For renormalization we performed runs at 
the same gauge couplings on 12 4 lattices up to (3 = 3.7, on 16 4 lattices up to 
= 4.05 and on 24 3 -32 at (3 = 4.1. We collected up to 2000 HMC trajectories 
for each finite temperature, and around 500 for each zero temperature point. 



3. Results 

The first quantity we study is the chiral condensate, ipip = (T/V)d/ dm log Z. 
This can be renormalized using the zero temperature condensate ipipo (this 
observable was studied in (32j|): 

m R ijjip R /ml = m(^V - V^o)/«4 ( 3 ) 

The renormalized condensate is plotted in Figure [2] together with our stag- 
gered estimate. The temperatures are converted to MeV again by using the 
physical value of Wq. We can see that there is still some lattice spacing de- 
pendence, but the N t = 8 results are very close to the staggered ones. One 
observes a broad cross-over, similar to the staggered results at physical quark 
masses 
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We have also determined the chiral susceptibility 

X ^ = (T/V)d 2 /dm 2 logZ, (4) 

but at the present level of our statistics renormalization resulted in large 
errors. Therefore we only show the bare susceptibilities in Figure [3J Even 
though the results obtained at our two lattice spacings cannot be compared 
directly, the peaks nicely signal the transition and we can see that the tran- 
sition temperature defined from this quantity has small lattice spacing de- 
pendence. 

The next quantity we study is the Polyakov loop. The bare quantity has a 
multiplicative divergence of the form exp[i 7 o(/3)/T] where the divergent term 



Fq can be determined up to a constant 34j. Different constants correspond 
to different renormalization schemes. We determine Fq entirely from finite 
temperature simulations in the following way. We perform runs with six 
different N t values, N t = 4, 5, 6, 7, 8 and 9, all lattice extents are 16 3 • N t . 
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Figure 2: The renormalized chiral condensate as a function of temperature on Nt = 6 
and 8 lattices. The upper temperature scale is for illustration only and it is based wq = 
0.1755 fm [3(|. The gray band shows our staggered estimate based on Nt = 6,8 and 10 
simulations. For details see |Appcndix B| 
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Figure 3: The bare chiral susceptibility as a function of temperature on our Nt — 6 and 
lattices. 
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Figure 4: The renormalized Polyakov loop as a function of temperature on N t = 6 and 
8 lattices as well as the staggered result. The black circle represents our renormalization 
condition, L R (T = 208MeV) = 1. 



We choose a fixed physical temperature such that these N t values span our f3 
range. This corresponds to a temperature of 208 MeV (having Nt = 9 at (3 = 
4.1). From these runs we can determine Fq(^) = 1/N t ■ logL at six /3 values. 
This can then be extended by interpolation to all of our couplings. This 
renormalization scheme corresponds to the condition Lr(T = 208MeV) = 1. 
The renormalized Polyakov loop is then given as: 

L R = L Q e- N ^\ (5) 

where L is the bare Polyakov loop. The result is shown in Figure |H We can 
see almost no lattice spacing dependence and an excellent agreement with 
the staggered results. 

Our final observable is the isospin susceptibility, 

Xi = {T/V)d 2 /d^\ogZ\^ =Q , (6) 

where /i/ is the isospin chemical potential, i.e. the quark chemical potentials 
are \i u = //// '2 and Hd = — //j/2. Obtaining results at non- vanishing chemi- 
cal potentials is very CPU demanding (see e.g. Even though our 



analysis to /x > is beyond the scope of the present Letter, including the 
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N t 


4 


6 


8 


10 


12 


£ = 2 overlap 


1.700 


1.588 


1.362 


1.241 


1.186 


£ = oo overlap 


1.619 


1.513 


1.290 


1.170 


1.117 


£ = 00 staggered 


2.235 


1.861 


1.473 


1.266 


1.164 


C, = 00 Wilson 


4.168 


2.258 


1.521 


1.265 


1.161 



Table 1: Stcfan-Boltzmann limits of the quark number susceptibility for three colors of 
overlap quarks with too = 1.3 for aspect ratios (£) of 2 and infinity. As a comparison we 
also give the infinite volume values for Wilson and staggered quarks. 



chemical potential even on the level of eq. (JS} is quite interesting. The rea- 
son is that there is an ongoing discussion in the literature about the proper 



inclusion of the chemical potential in the overlap operator [37], |38|, |39(. We 
follow Reference 37J and define the chemical potential as a fourth, imaginary 
component of the temporal gauge field and use the generalization of the sign 
function: sgn(z) = sgnRe(2i). The second derivative can be calculated using 
the formulas of Appendix A As a tree level improvement we normalize all 
susceptibilities with the corresponding Stefan-Boltzmann (SB) values which, 
for our choice of mo = 1.3 are given in Table [1] for a number of NtS both for 
infinite volume and our relatively small aspect ratio, £ = 2. Our results for 
the quark number susceptibility are shown in Figure [5j Again, the N t = 8 
results are very close to the staggered ones. 



4. Conclusions, outlook 

We presented results for the temperature dependence of several observ- 
ables using dynamical overlap fermions. Our results show that on Nt — 6 and 
8 lattices cutoff effects are still present but not severe. The comparison with 
staggered results obtained on N t = 6,8 and 10 lattices is quite encouraging: 
there is a good chance for a reliable continuum extrapolation from Nt 10 
lattices with overlap fermions. 

There are several remaining issues to investigate in future studies. A 
particularly interesting question is the dependence of the results on the global 
topology and how it disappears in the thermodynamic limit. Such analysis 
requires a series of runs on larger volumes. To approach the continuum limit 
one needs at least one more lattice spacing. Finally, one might also include 
the strange quark and decrease the light quark masses to their physical values. 
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Figure 5: The isospin susceptibility as a function of temperature on N t = 6 and 8 lat- 
tices. The gray band again shows the corresponding staggered result. All data sets were 
normalized by their corresponding SB limits (cf. Table [1]). 
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Appendix A. Derivatives of the Zolotarev approximation 

Here we present results on the first and second derivatives of the sign 
function. In the following the S symbol stands for a general derivation oper- 
ator, for example a derivation with respect to the gauge links or with respect 
to the chemical potential. The Zolotarev approximation of the sign function 
of an operator h can be written in a partial fraction expansion as: 

sgn(/i) « h I c + T^J— ) = h ( c o + J2 CiQi ) ' (A- 1 ) 
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where Qi = (h 2 + 1 was introduced. The first derivative is 

8sgn(h) « <J/i ^c + ^ QQij - CiQih5h 2 Qi, (A.2) 



whereas the second derivative is 



5 2 sgn(/i) ^5 2 h[c + J2 ) + J2 Qi [25h 2 Q t h5h 2 - h5 2 h 2 - 25h5h 2 ] c t Qi 

(A.3) 



Appendix B. Details of the staggered calculations 

In this appendix we describe our staggered analysis. The sole goal of this 
study was to provide a basis of comparison for the overlap data. To suppress 
most of the taste-breaking effects of the staggered action we used four levels 
of stout smearings (p = 0.125) in the fermionic sector, which contained here 
the two light quarks only. We determined a two-flavor LCP for this staggered 
action using the same definition as in the case of the overlap simulations. 

To this end we simulated sixteen ensembles (four lattice spacings with 
four quark masses each). The smallest mass was always close to the final 

LCP's value. In the range (3 = 3.8 4.1 we could fit the emerging am q (f3) 

and wo/a(P) functions to sub-percent accuracy. 

For the renormalization of the chiral condensate we made a chiral interpo- 
lation (or extrapolation for some lattice spacings) of the vacuum condensate 
(#)o- 

The Polyakov loop requires renormalization, too, here we calculated the 
static potential for our sixteen T = ensembles. One can select any physical 
distance to remove all divergences. Since wo is our most accurately known 
scale, we used \/8w which is around 0.5 fm and defined Z = exp(V(r = 
y/8w )/2). After this renormalization our final estimate for the Polyakov loop 
was multiplied by an exp(— S/(wqT)) function, which facilitates finite trans- 
formations between different renormalization schemes. We set 5 = 0.1566 
to approximately reproduce the renormalization condition used in the over- 
lap results. The renormalization procedure might seem elaborate, yet the 
dominant source of error comes from our finite temperature statistics. 

Our finite temperature simulations were performed with the same aspect 
ratio as with the overlap action, but here, in addition to N t = 6 and 8, we 
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also made a set of N t = 10 ensembles. In most of the temperature range we 
had two lattice spacings. The agreement between the results at these lattice 
spacings was even slightly better than in our previous work in full QCD 34 . 
For our "staggered estimate" we use the Nt = 10 at low temperatures and 
we present N t = 8 data starting from a temperature where both N t = 8 and 
10 data were present and they were in agreement. 
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